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Abstract 

The  modeling  of  proton  exchange  membrane  fuel  cells  (PEMFC)  may  work  as  a  powerful  tool  in  the  development  and  widespread  testing  of 
alternative  energy  sources  in  the  next  decade.  In  order  to  obtain  a  suitable  PEMFC  model,  which  can  be  used  in  the  analysis  of  fuel  cell-based 
power  generation  systems,  it  is  necessary  to  define  the  values  of  a  specific  group  of  modeling  parameters.  In  this  paper,  the  authors  propose  a 
dynamic  model  of  PEMFC,  the  originality  of  which  lays  on  the  use  of  non-integer  derivatives  to  model  diffusion  phenomena.  This  model  has  the 
advantage  of  having  least  number  of  parameters  while  being  valid  on  a  wide  frequency  range  and  allows  simulating  an  accurate  dynamic  response 
of  the  PEMFC. 

In  this  model,  the  fuel  cell  is  represented  by  an  equivalent  circuit,  whose  components  are  identified  with  the  experimental  technique  of  electro¬ 
chemical  impedance  spectroscopy  (EIS).  This  identification  process  is  applied  to  a  commercially  available  air-breathing  PEMFC  and  its  relevance 
is  validated  by  comparing  model  simulations  and  laboratory  experiments.  Finally,  the  dynamic  response  derived  from  this  fractional  model  is 
studied  and  validated  experimentally. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

It  is  believed  that  there  will  be  a  time  in  the  future  when  global 
energy  demands  will  be  met  by  some  sources  of  energy  other 
than  fossil  fuels.  Thus,  fuel  cells,  in  particular  proton  exchange 
membrane  fuel  cells  (PEMFC),  are  expected  to  play  a  major  role 
in  the  future  energy  sector.  PEMFC  are  particularly  attractive 
for  use  in  vehicles  as  a  replacement  to  the  internal  combustion 
engines.  They  also  seem  to  be  a  promising  source  to  be  used 
in  residences,  industries  and  small-  and  large-scale  distributed 
generation  systems.  The  low  operating  temperature  of  a  PEMFC 
(typically  <90  °C)  allows  easy  start-up  and  fast  response  to  load 
variations  and  operating  conditions.  Nevertheless,  several  issues 
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need  to  be  resolved  before  fuel  cells  can  be  commercially  viable. 
Indeed,  the  need  of  precise  water  management,  the  dehydration 
of  membrane,  the  complex  electrode  kinetics,  the  mass  transport 
and  the  slow  rate  of  oxygen  reduction  are  the  most  significant 
limiting  factors  on  the  fuel  cell  performances. 

In  order  to  better  understand  the  physical  phenomena  in  the 
fuel  cells  and  describe  their  steady  state  and  dynamic  behaviors, 
the  modeling  of  fuel  cells  has  become  more  and  more  important 
over  the  years  in  order  to  simulate  the  behaviour  of  the  fuel  cell 
integrated  in  a  power  system.  Such  model  must  be  sufficiently 
complex  to  take  into  account  all  the  electrochemical  phenomena 
and,  at  the  same  time,  should  be  able  to  be  integrated  into  a 
complete  system  simulation.  Furthermore,  these  models  must  be 
modular  so  as  to  allow  easily  the  testing  of  various  technological 
solutions. 

In  this  article,  the  authors  have  focused  on  the  dynamic  mod¬ 
eling  of  a  fuel  cell.  Such  modeling  to  simulate  the  transient 
response  of  a  PEMFC  is  studied  only  recently.  In  the  past,  the 
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Nomenclature 

Nomenclature 

b  Tafel  slope  (V  dec-1) 

Cdi  double  layer  capacitance  (F  m-2) 

C  concentration  of  species  (molm-3) 

D  diffusion  coefficient  (m2  s-1) 

E  potential  (V) 

F  Faraday’s  constant  (96,500  C  mol-1) 

j  current  density  (A  m-2) 

jo  exchange  current  density  (A  m-2) 

j\  limiting  current  density  (A  m-2) 

nQ  total  number  of  exchanged  electrons 
P  anode  and  cathode  gas  pressure  (Pa) 

R  ideal  gas  constant  (8.314  J  mol-1  K-1) 

Rm  membrane  resistance  (Q  m2) 

s  Laplace  operator  =  i  x  co 

T  temperature  (K) 

fA cell  cell  potential  (V) 

Z  impedance  (£2) 

Greek  symbols 

8  gas  diffusion  layer  thickness  (m) 

£  porosity 

y  roughness  factor  (m2  m-2) 

0  overpotential  (V) 

co  angular  frequency  (rad  s - 1 ) 

r  time  constant  (s) 

Subscripts  and  superscripts 
a  anode 

c  cathode 

k  either  anode  or  cathode 

m  membrane 

*  referring  to  the  gas  channel/gas  diffusion  layer 

interface 

eff  referring  to  the  effective  (apparent)  value 


mathematical  models  of  the  PEMFC  have  usually  been  restricted 
to  the  steady  state  conditions.  Bernardi  [1]  first  proposed  a  one¬ 
dimensional  model  in  order  to  study  the  water  management  and 
to  identify  the  humidification  conditions  which  induce  either 
the  dehydration  of  the  membrane  or  excessive  flooding.  The 
sensitivity  to  the  water  balance  of  the  PEMFC  performance 
was  demonstrated.  Some  other  models  were  derived  from  the 
original  one  to  take  into  account  heat  management  [2],  mass 
transport  in  the  gas  diffusion  electrode  [3]  or  to  introduce  a  dif¬ 
ferent  treatment  of  the  electrochemical  reaction  [4].  Bernardi  et 
al.  [5]  and  Springer  et  al.  [6]  also  presented  a  one-dimensional 
model  to  investigate  the  factors  that  limit  cell  performance  and 
to  elucidate  the  mass  transport  mechanism  within  the  complex 
network  of  gas,  liquid  and  solid  phases  constituting  the  gas  dif¬ 
fusion  electrode.  Contrary  to  the  previous  models,  this  latter 
one  considered  the  molar  changes  along  the  gas  channel.  These 
approaches  made  it  possible  to  evaluate  the  losses  in  the  cathode 
gas  backing  and  the  catalyst  layer. 


So  far  efforts  have  been  made  to  model  the  problem  in 
three  dimensions.  A  higher  dimensionality  makes  it  possible 
to  describe  the  hydrodynamics  and  multi  component  transport 
inside  the  flow  channel  for  reactant  distribution.  Multi  com¬ 
ponent  models  were  thus  developed  for  cathode  or  for  whole 
PEMFC,  for  both  single-phase  [7]  and  two  phase  flow  [8].  These 
models  showed  that  the  gas  distributor  geometry  has  a  signif¬ 
icant  effect  on  the  diffusion  of  the  reactants  and  products  in 
determining  the  performance  of  the  cell.  Moreover,  they  con¬ 
firmed  that  the  performance  of  cathode  is  strongly  influenced  by 
the  presence  of  liquid  water  and  its  removal  rate,  especially  at 
high  current  density. 

The  dynamic  study  of  a  fuel  cell  appears  to  be  of  great  interest 
to  provide  detailed  understanding  of  mass  and  charge  transport 
through  the  gas  diffusion  electrodes  and  is  of  extreme  impor¬ 
tance  for  the  control  strategy  and  system  management  in  power 
generation  systems,  especially  when  there  are  power  injections 
into  the  network.  Thus  a  first  equivalent  circuit  was  proposed  to 
simulate  the  impedance  spectrum  of  a  PEMFC  [9] .  This  study 
showed  that  the  impedance  of  a  fuel  cell  is  a  powerful  tool  in 
order  to  characterize  the  intra-electrode  processes  occurring  in 
gas  diffusion  electrodes. 

Some  dynamic  models  of  PEMFC  [10-13]  have  been  devel¬ 
oped  based  on  the  physical  and  chemical  knowledge  of  the 
phenomena  occurring  inside  the  cell.  These  models  are  gen¬ 
erally  implemented  in  Simulink-Matlab  environment.  Results 
indicated  that  the  transient  response  of  the  PEMFC  to  reach 
steady  state  is  less  than  10  s.  An  important  effect  of  water  man¬ 
agement  was  exhibited  too  [12]. 

Some  other  models  were  developed  in  a  more  “System” 
approach  [14].  These  models  are  relatively  simple  which  allow 
correctly  simulating  the  behavior  of  the  fuel  cell  inserted  in  an 
electrical  network.  The  main  disadvantage  of  these  models  is 
that  they  are  far  away  from  the  physics  of  the  fuel  cell,  i.e. 
their  parameters  are  physically  non-representative.  In  parallel, 
some  other  multi-scale  models  were  elaborated  [15]  to  predict 
the  dynamic  and  steady  state  behaviors  for  the  triple  contact 
(base  of  the  electrochemical  reaction).  These  models  made  it 
possible  to  describe  quantitatively  the  reactional  mechanisms, 
the  polarization  curves  and  impedance  spectra  of  the  fuel  cell. 

In  this  work,  the  authors  have  linked  these  two  approaches 
(components  and  systems)  by  using  the  fractional  approach  in 
order  to  obtain  an  adequate  model  of  fuel  cell  impedance,  which 
can  be  used  in  a  system  simulation.  Non-integer  derivation  has 
already  been  used  to  correctly  model  the  diffusion  phenomenon 
of  magnetic  field  in  electrical  machines  [16]  which  is  a  phys¬ 
ical  phenomenon  similar  to  one  found  in  the  electrochemical 
devices.  The  resulting  models  are  precise,  having  less  num¬ 
ber  of  parameters  and  being  valid  on  a  wide  frequency  range. 
Moreover,  parameters  of  such  non-integer  order  models  have 
a  close  link  with  the  physical  characteristics  of  the  machines. 
This  method  is  very  useful  for  the  optimal  real-time  control  of 
PEMFC  running  on  a  load. 

In  this  study,  the  authors  have  thus  used  non-integer  deriva¬ 
tives  to  model  the  diffusion  phenomenon  of  gases  on  electrodes. 
Then,  this  new  modeling  technique  is  compared  with  the  clas¬ 
sical  one  while  taking  into  account  the  characteristic  response 
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of  the  fuel  cell  (i.e.  polarization  curves  and  impedance  spec¬ 
tra).  In  order  to  identify  the  model  parameters,  a  complete  set 
of  experimental  results  was  obtained  for  a  commercially  avail¬ 
able  air-breathing  PEMFC  using  the  technique  of  impedance 
spectroscopy.  Finally,  the  dynamic  response  derived  from  this 
fractional  model  is  studied  and  validated  experimentally. 

2.  Modeling  of  fuel  cell 

2.7.  Hypotheses 

The  air-breathing  PEMFC  system  is  considered  as  isothermal 
and  isobaric.  Both  these  approximations  appear  to  be  valid  since 
these  conditions  are  normally  achieved  in  a  small  single-cell 
experimental  investigation.  The  total  pressure  at  the  anode  and 
cathode  compartment  and  within  the  gas  diffusion  electrodes  is 
considered  as  constant.  Moreover,  ionic  ohmic  drop  in  the  active 
layer  and  the  electronic  ohmic  drop  in  the  current  collectors  can 
be  neglected  owing  to  the  high  electronic  and  ionic  conductiv¬ 
ities,  and  thus  lead  to  the  absence  of  voltage  drop.  It  is  also 
supposed  that  the  gas  permeation  in  the  membrane  is  negligible. 

2.2.  Modeling  in  steady- state  regime 

The  physical  model  used  to  obtain  the  following  governing 
equations  for  steady  state  and  transient  responses  of  PEMFC 
is  based  on  a  previous  d.c.  and  a.c.  model  [15].  Fig.  1  shows  a 
schematic  illustration  of  a  self-breathing  PEMFC  that  consists  in 
a  membrane  sandwiched  between  two  gas  diffusion  electrodes, 
this  assembly  being  pressed  between  two  current  collectors  and 
end  plate.  A  simplified  view  of  an  electrode  and  the  induced 
variation  of  the  concentration  of  reactant  gazes  along  the  width 
of  diffusion  layer  are  also  presented. 

In  Fig.  1,  8  represents  the  diffusion  layer  thickness,  C|  = 
Pk/RT  is  the  concentration  of  oxygen  (O2)  or  hydrogen  (H2) 
in  the  gas  channels,  C*  the  concentration  of  hydrogen  or  oxy¬ 


gen  at  the  active  layer/diffusion  layer  interface  and  7\  the  gas 
pressure. 

Far  from  the  equilibrium,  the  hydrogen  oxidation  and  oxygen 
reduction  are  classically  described  by  the  laws  of  Tafel  [15]: 


Ja  =  JOa  exp 


and 


jc  —  JOc  exp 


where  jok  is  the  exchange  current  density,  bk  the  Tafel  slope, 
1*7*1  =  |  Ek  —  I  absolute  value  of  the  overpotential,  R  the 
ideal  gas  constant  and  T  is  the  temperature.  The  subscript 
represents  either  anode  (a)  or  cathode  (c). 

It  has  been  shown  in  [15]  that  the  ratio  of  the  concentration 
of  the  gas  at  the  interface  is  closely  linked  to  the  current  density 
as  follows: 


Ck  _  1  jk 

cf~  ~Ji 


Finally,  the  equations  (1)  are  expressed  using  the  limiting  current 
densities  and  introducing  the  roughness  factor  of  the  electrodes: 


and 

(3) 

While  the  limiting  current  density  can  be  defined  as: 

Deff 

M  =  -f-C*kne,kF  (4) 

where  y*  is  the  roughness  factor,  nQ k  the  number  of  exchanged 
electrons,  F  the  Faraday’s  constant,  Dff  =  sl -5 Dk  the  effective 
gas  diffusion  coefficient  and  s,  the  porosity  of  electrodes.  The 
coefficient  y*  is  introduced  to  account  for  the  roughness  property 


J a  —  Pa  JOa  exp 


Jc  —  Yc  J Oc  exP 


2.31/7 


a  I 


^a 

2.3 1 77c  | 


b( 


1  - 


1  - 


Ja 

Jla 

Jc 

Jlc 


h2 


cathode  Electrode  compartment 

Assembly 


Fig.  1.  Simplified  diagram  of  an  electrode  of  PEMFC. 
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of  the  porous  active  layer.  It  means  that  the  catalyst  particles  are 
distributed  in  all  the  volume  of  the  active  layer. 

From  Eq.  (2),  it  is  then  useful  to  express  each  electrode  poten¬ 
tial  in  terms  of  current  density: 


Ek  —  E\\yk  + 


jk 


Yk  •  jok 


where  E^k  is  the  thermodynamic  potential  of  the  electrode.  Cal¬ 
culating  the  overpotentials  of  cathode  and  anode  and  fixing  the 
values  of  the  electrochemical  parameters,  the  polarization  curve, 
i.e.,  Uce ii  as  a  function  of  current  density  j,  can  then  be  plotted. 
Thus,  by  subtracting  the  ohmic  drop  ( Rm-j )  and  the  absolute  val¬ 
ues  of  overpotentials  of  cathode  and  anode  from  the  reversible 
thermodynamic  potential  of  the  fuel  cell  (1.23  V),  the  polariza¬ 
tion  curve  is  given  by 


f^cell  —  1.23  \t]sl\  I  Oc  I  ’  j  (6) 

2.3.  Modeling  in  dynamic  regime 


then  given  by: 

^Total(^)  =  Za(<X>)  +  Rm  +  ZC(<X>)  (8) 


Around  a  stationary  operating  point,  we  obtain  the  charge  trans¬ 
fer  resistance  R±  by  differentiating  the  Eq.  (3)  with  respect  to 
the  over  potential  (77): 


1 

Rtk  =  - 

djk/dOk 


1 


Ykjok  |f  exp 


(  2.3|>»-|  \ 

l  bt  ) 


while  the  analytical  impedance  of  mass  diffusion  (i.e.  Warburg 
impedance)  is  expressed  as  [15,17]: 


Zwk(s)  = 


djk/dCk  8  tanh(0Tp 
djk/dm  neFDf  ^k 


(10) 


where  s  =  I,a>  is  the  Laplace  operator,  and  r k  =  S2/Df,  the  time 
constant  of  diffusion. 

Here,  Eq.  (10)  can  be  simplified  introducing  the  non-linear 
term  Ak(j)  which  depends  on  the  current 


For  the  dynamic  and  small  signal  mode,  the  same  assumptions 
have  been  conceived  as  those  supposed  in  the  stationary  regime. 


2.3.1.  Analytic  electrical  model  of  fuel  cell 

When  there  is  no  mass  transport  limitation,  the  redox  reaction 
is  simply  represented  by  an  equivalent  electrical  circuit  of  paral¬ 
lel  RC  cells.  However,  when  there  are  considerable  variations  of 
the  interfacial  concentrations  on  electrodes,  the  redox  reaction 
is  represented  by  an  equivalent  circuit  of  parallel  ZfC  [15].  The 
faradic  impedance  Zf  is  composed  of  two  impedances  in  series: 
a  charge  transfer  resistance  Rt  and  an  impedance  of  diffusion  of 
the  reduced  species  on  cathode  side  or  oxidized  species  on  anode 
side.  This  impedance  of  diffusion  is  called  Warburg  impedance 
and  is  represented  by  Zw  [15,17]. 

The  impedance  of  an  electrode  then  corresponds  to  the  paral¬ 
lel  combination  of  the  faradic  impedance  Zf  and  the  double  layer 
capacitance  (Cdi)  to  account  for  the  dynamics  of  the  changing 
concentration  in  the  gas  backing  layer  and  the  charge  stored 
in  the  interfacial  capacitance.  The  expression  of  the  electrode 
impedance  is  thus  given  by: 


Z(co)  = 


Zfjco) 

1  +  iooC&xZfoo) 


Finally,  the  total  impedance  of  the  fuel  cell  is  composed  of 
two  impedances,  one  impedance  for  each  electrode  (anode  and 
cathode),  in  series  with  the  internal  resistance  Rm  linked  to  the 
membrane  (Fig.  2).  The  total  impedance  Zpotai  of  the  fuel  cell  is 


Zwk(s)  =  Ak(j) 


tanh  (y/sn) 

V^k 


where  Ak(j) 


1  5 

I" 


(ii) 


The  effective  double  layer  capacitance  C^f  is  defined  as  C^[f  = 
y  •  Cdu  where  y  represents  the  roughness  factor  whose  value  is 
of  the  order  of  100. 


2.3.2.  Distributed  parameter  modeling  approach 

The  presence  of  the  tangent  hyperbolic  function  in  the  ana¬ 
lytical  expression  of  the  Warburg  impedance  does  not  enable  us 
to  draw  an  equivalent  electrical  circuit  which  is  directly  usable 
to  model  the  transient  response  of  the  fuel  cell  in  any  operating 
mode. 

The  classical  approach  of  modeling  Warburg  impedance  con¬ 
sists  in  the  decomposition  of  the  tangent  hyperbolic  function  in 
the  form  of  a  series  (theoretically  infinite)  and  in  the  identi¬ 
fication  of  these  serial  terms  of  the  RC  cells  parameters.  The 
principal  drawback  of  this  approach  is  that  it  required  a  trun¬ 
cation  if  the  infinite  series.  The  order  of  truncation  is  chosen 
according  to  a  certain  tolerable  error.  Currently,  it  is  the  major 
problem  related  to  the  modeling  of  systems  with  distributed 
parameters  of  infinite  order. 

The  decomposition  of  the  function  tanh  in  the  form  of  a  series 
(series  of  Foster)  is  given  by  [16]  (Fig.  3): 

tanh(v)  =  2x 


The  Warburg  impedance  (Eq.  (11))  then  becomes: 


00 

E— 

n=\ X2  + 


1 


n(2n—  1) 


i2 


00 

Zwk(s)  =  2  Ak(j)^2 - 

n=lSTk  + 


1 


n(2n  —  1) 
2 


(12) 


Fig.  2.  Complete  equivalent  circuit  of  the  fuel  cell. 
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R,  R2 


r#-i 

w — 

C 

/ 

C 

7 

c, 

Fig.  3.  Foster  series  equivalent  of  Warburg  impedance. 


We  can  then  model  this  Warburg  impedance  by  infinite  number 
of  parallel  RC  cells  connected  in  series: 

The  expression  of  resistance  Rn  and  capacitance  Cn  for  each 
cell  n  is  given  by 


2  A(j) 


n(2n  —  1) 
2 


and  Cn 


-r 

n(2n  —  \) 

2 

l 

2 

2  A(j) 

It  should  be  noted  that  these  parameters  depend  on  the  current 
density  j. 


2.3.3.  Non-integer  order  modeling 

Non-integer  derivatives  are  introduced  in  Eq.  (11)  using 
another  approximation,  i.e.,  the  second-order  Taylor’s  expan¬ 
sion  of  the  tanh  function  [16]: 

tanh  (ffsr)  sinh(,v/vr)  1  \fsr  1  1 

- , -  -  - , - ,  - ,  , 

ff~SX  cosh(^sr)  ffsX  1  +  sdf  ffsX  y/\  +  sx 

(13) 

This  limited  expansion  however  remains  valid  at  high  frequen¬ 
cies;  indeed,  using  the  approximation  of  tanh  function  for  high 
frequencies: 

tanh(V^r)  sinh(Vvr)  1  e^  11  1 

- , -  -  - _  —  - ,  ^  , 

ff~sx  cosh(^sr)  ^srooev^r  ffsxoo  ffsroo  _j_  sr 

(14) 

By  using  the  approximation  (13)  we  can  simplify  the  analytical 
expression  of  the  Warburg  impedance  given  by  the  Eq.  (11),  as 
follows: 

AkU) 

Zwk(s)  =  .  (15) 

V  1  +  srk 


The  half-order  fractional  model  of  the  Faradic  impedance  is  then 
given  by 

Ak(j ) 

Zfk(s)  =  Rtk  +  Zwk(s)  =  Rtk  +  (16) 

Vl  +STk 

The  total  impedance  of  the  electrode  is  then  given  by 


^electrode  C?)  — 


1 


Zfk(s) 


-L 


(17) 


Finally,  the  total  impedance  of  the  fuel  cell  is  as  usually  given 
by  the  relation  (8). 


2.4.  Comparison  of  different  models 

In  Fig.  4,  the  various  models  of  Warburg  impedance,  pre¬ 
sented  in  Section  2.3,  are  compared  in  frequency  domain,  i.e., 

•  Analytical  Warburg  impedance  given  by  Eq.  (11). 

•  Classical  model  given  by  Eq.  (12)  using  20  RC-cells  in  series. 

•  Half-order  fractional  model  given  by  Eq.  (14). 

These  curves  have  been  simulated  for  A(j)  =  1 .0  £2  rad1/2  and 
r  =  82/D  =  1  s.  Thus,  it  can  be  noted  that  the  difference  between 
the  fractional  Warburg  impedance  and  the  analytical  Warburg 
impedance  is  very  small  over  a  wide  range  of  frequencies.  This 
small  difference  that  exists  between  the  two  models  can  be 
related  to  the  second-order  expansion  of  the  tangent  hyperbolic 
function.  This  limited  expansion  however  remains  valid  at  high 
frequencies. 

It  is  also  worth  mentioning  that  the  fractional  Warburg 
impedance  represents  resistive  behavior  at  low  frequency,  which 
is  easily  observable  by  a  horizontal  line  (Fig.  4a)  and  is  charac¬ 
terized  by  a  zero  phase  (Fig.  4b).  Moreover,  asymptotically,  the 
Bode  plot  of  fractional  model  is  a  straight  line  having  a  slope 
of  —  lOdB  per  decade;  while  its  phase  is  constant  and  equal  to 
—45°.  The  asymptotic  behavior  of  Warburg  impedance  is  then 
totally  taken  into  account  by  the  fractional  model  contrary  to  the 
classical  one,  which  behaves  like  a  capacitance  for  high  frequen¬ 
cies.  The  number  of  parameters  describing  each  model  can  also 
be  highlighted:  2  for  fractional  model  and  40  for  the  classical 
R-C  one! 


Fig.  4.  Comparison  of  Warburg  impedance,  analytical  expression  (••••),  classical  model  using  20  RC  cells  ( — ),  fractional  model  ( — ).  (a)  Variation  of  gain  of 
impedance,  (b)  Variation  of  phase. 
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Fig.  5.  Comparison  between  different  models  of  impedance  (in  Ohms)  for  j  =  1 .2  A  cm  2 ,  analytical  expression  ( — ),  classical  model  using  20  RC  cells  (••••),  Fractional 
model  ( — ).  (a)  Total  impedance  of  the  PEMFC.  (b)  Zoom  on  the  diffusion  arcs. 


Similarly,  the  total  impedance  of  the  fuel  cell,  using  the 
two  different  models  of  Warburg  impedance  can  be  compared 
(Fig.  5). Therefore,  Fig.  5  shows  that  the  classical  model  with 
(R-C)  cells  seems  to  be  more  precise  than  the  fractional  model, 
particularly  in  the  range  of  frequencies  where  the  diffusion  phe¬ 
nomenon  appears,  i.e.,  0.1  and  100  Hz  (Figs.  4  and  5).  However, 
the  overall  shape  of  the  fractional  model  is  in  a  good  agree¬ 
ment  with  the  shape  obtained  from  the  analytical  impedance. 
Although  there  exists  a  mean  error  of  0.18%  compared  to  the 
analytical  expression  on  a  given  frequency  range,  this  error  is 
negligible  for  the  reduced  number  of  parameters  of  this  model. 
However,  it  is  important  to  note  that  this  error  depends  on  the 
current  density  and  that  its  value  increases  with  current  density 
but  remains  lower  than  0.5%. 

In  spite  of  this  inconvenience,  the  important  point  to  be  noted 
is  that  the  number  of  parameters  of  the  fractional  model  is 
reduced.  Indeed,  the  classical  (R-C)  model  has  2  x  N  param¬ 
eters  for  N  cells  (R-C)  in  series,  while  the  fractional  model  is 
characterized  only  by  2  parameters.  It  is  a  consequence  of  the 
intrinsic  compactness  of  the  fractional  systems.  This  property, 
which  leads  to  system  order  reduction,  allows  decreasing  the 
simulation  times  of  PEMFC  system. 

3.  Validation  of  non-integer  order  model 

3.1.  Experimentation 

Impedance  measurements  were  carried  out  on  a  commer¬ 
cially  available  air-breathing  PEM  fuel  cell  provided  by  PAX- 
ITECH.  It  contained  a  Nafion®-115  membrane.  The  platinum 
loading  on  the  electrodes  was  0.3  mg  cm-2.  The  electrode  sur¬ 
face  area  was  25  cm2  and  the  diffusion  layer  thickness  of 
200  juim.  The  experimental  results  were  obtained  at  ambient  tem¬ 
perature  (30  °C)  while  the  cell  functioned  with  air  at  cathode  and 
hydrogen  at  anode.  The  test  bench  consisted  of  a  Solartron  1250 
frequency  response  analyzer  coupled  to  a  Solartron  1286  elec¬ 
trochemical  interface.  EIS  experiments  were  performed  under 
galvanostatic  mode  with  modulating  AC  current  amplitude  of 


100mA  and  range  of  the  frequency  lies  between  100  mHz  and 
65  kHz  with  10  points  per  decade. 

In  order  to  study  the  fuel  cell,  impedance  measurements  were 
made  at  open  circuit  and  at  various  points  along  the  polarization 
curve  by  changing  the  load  resistance  in  order  to  vary  the  output 
current  of  the  cell  [18,19].  All  the  impedance  spectra  shown 
in  Fig.  6  are  modified  with  respect  to  the  corresponding  load 
resistance. 

From  the  impedance  diagrams  of  fuel  cell,  it  is  possible  to 
determine  a  set  of  parameters.  The  high  frequency  intercept  at 
real  axis  corresponds  to  the  value  of  the  internal  resistance  R[  of 
which  a  fraction  corresponds  to  the  membrane  resistance,  Rm ; 
here,  it  varies  approximately  between  0.23  and  0.32  Q.  More¬ 
over,  the  low  frequency  limit  is  the  polarization  resistance  Rp, 
which  is  defined  as  the  derivative  (i.e.  the  slope)  of  the  polariza¬ 
tion  curve  of  the  cell.  In  Fig.  6,  from  diagram  a  to  diagram  k, 
Rp  varies  approximately  from  3  to  0.4  Q,  and  then  it  becomes 
almost  constant  at  0.4  Q. 

In  Fig.  6b  and  c,  the  diagrams  of  impedance  are  composed 
of  only  single  arc.  This  shows  that  at  low  current  densities, 
the  impedance  of  fuel  cell  is  dominated  by  the  phenomenon 
of  charge  transfer.  At  low  frequencies  (i.e.  for  the  frequencies 
lower  than  0.5  Hz)  the  diagrams  of  impedance  do  not  have  a 
particular  form,  and  we  cannot  give  their  physical  significance. 
In  this  study,  we  will  ignore  these  points  corresponding  to  the 
frequencies  lower  than  0.5  Hz.  Moreover,  it  is  well  known  that 
for  high  current  densities,  the  diagrams  of  impedance  are  classi¬ 
cally  made  of  two  arcs:  the  first  one,  the  high  frequency  arc,  is  an 
arc  of  charge  transfer;  while  the  second  one  is  an  arc  related  to 
mass  transport  phenomenon  in  the  diffusion  layer  (mainly  due 
to  the  diffusion  of  oxygen  at  cathode).  Moreover,  it  can  be  stated 
from  literature  [20]  that  the  arc  related  to  charge  transfer  (i.e.  the 
high  frequency  arc)  decreases,  whereas  the  arc  related  to  mass 
diffusion  (i.e.  the  low  frequency  arc)  increases,  with  the  current 
density.  However,  in  our  measurements,  these  two  arcs  cannot 
be  observed  separately.  This  can  be  explained  on  the  one  hand 
by  the  fact  that  the  output  current  of  the  cell  is  not  sufficiently 
high  and  on  the  other  hand  by  a  high  value  of  the  double  layer 
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Fig.  6.  Impedance  diagrams  of  the  fuel  cell  measured  along  the  polarization  curve,  (a)  Open  circuit,  (b)  /=  0.564  A  and  V  =  0.584  V.  (c)  1=  0.906  A  and  V  =  0.480  V. 
(d)  1=  1.116 A  and  V=0.411  V.  (e)  7=1.28A  and  V=0.357V.  (f)  7=1.38A  and  V=0.315V.  (g)  7=1.512A  and  V=0.271  V.  (h)  7=1.687  A  and  V=0.244V.  (i) 
7=  1.752  A  and  V=  0.237  V.  (j)  7=  1.80  A  and  V=  0.233  V.  (k)  1=  1.84  A  and  V=  0.220  V. 


capacitance,  which  is  hiding  the  two  different  phenomena  into 
one  arc. 

3.2.  Identification  of  model  parameters 

The  non-linear  least  square  fitting  algorithm  of  MATLAB® 
was  used  to  identify  the  relevant  parameters  of  the  fractional 
model  of  PEMFC  developed  in  Section  2.3.  The  principle 
of  this  algorithm  is  to  minimize  the  square  of  a  nonlinear 
function  while  finding  the  best  values  of  the  unknown  vari¬ 
ables  (i.e.  model  parameters)  starting  from  their  given  initial 
values. 

It  should  be  noted  that  the  impedance  points  measured  at 
low  frequencies,  whose  physical  significance  is  not  clear,  and 
the  points  of  impedance  having  a  positive  imaginary  part  (due 


to  wiring  inductance)  are  not  taken  into  account  for  this  iden¬ 
tification  process.  As  shown  in  Section  2.3,  the  parameters  of 
the  fractional  model,  which  are  to  be  identified,  are  Rtk ,  A&,  r, 
Rm  and  y-Cdi-  The  effective  double  layer  capacitance,  y-Cdi  is 
determined  by  an  algorithm  of  identification  and  then  this  value 
is  used  to  do  “fitting”  for  the  other  parameters  of  the  model 
with  experimental  results.  The  value  of  the  time  constant  r  is 
calculated  from  the  frequency  corresponding  to  the  peak  of  the 
diffusion  arc.  As  far  as  internal  resistance  R[  is  concerned,  it  cor¬ 
responds  to  the  distance  between  the  origin  and  the  intersection 
of  the  impedance  spectrum  with  the  real  axis. 

The  remaining  parameters  (Rtk  and  Af)  are  identified  by  using 
non-linear  least  square  algorithm.  Table  1  shows  the  values 
of  the  model  parameters  identified  for  four  impedance  spectra 
corresponding  to  current  densities  ranging  between  0.068  and 
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Table  1 


Values  of  model  parameters  identified  for  different  current  densities 


Parameters 

j  =  0.068  A  cm  2 

j  =  0.070  A  cm  2 

7  =  0.072  A  cm  2 

j  =  0.074  A  cm  2 

Units 

Aa 

6.67  x  10“9 

3.48  x  10“9 

2.48  x  10“9 

5.53  x  10“10 

Q 

Ac 

1.80  x  10“5 

1.28  x  10“5 

1.10  x  10“5 

7.0  x  10"6 

Q 

Aa 

3.56  x  10-10 

1.23  x  10“8 

2.0  x  10“8 

2.62  x  10“8 

£2  rad 1/2 

Ac 

3.0  x  10"5 

4.6  x  10“5 

6.37  x  10“5 

8.4  x  10“5 

£2  rad 1/2 

A 

0.236 

0.254 

0.283 

0.313 

£2 

y-Cd  i 

152 

115 

115 

115 

F 

T 

0.0404 

0.0404 

0.0404 

0.0404 

s 

Mean  error 

1.41 

2.56 

5.50 

8.22 

% 

0.074  A  cm-2,  and  the  mean  error  corresponding  to  each  set  of 
parameters. 

As  expected,  the  four  parameters,  Ata,  Rtc ,  Aa  and  Ac,  turn 
out  to  be  more  sensitive  to  the  current  density  than  the  other 
ones  (r,  R[  and  y-Cdi).  The  evolution  of  the  model  parameters 
with  current  densities  shows  that  the  higher  the  output  current 
is,  more  the  fuel  cell  is  limited  by  the  diffusion  phenomenon  of 
the  gases  (increase  in  A&);  on  the  other  hand,  the  charge-transfer 
resistance,  Rtk ,  decreases  with  the  current  density.  Moreover,  it 
is  clearly  shown  that  the  limitation  by  the  diffusion  phenomenon 
at  cathode  is  more  significant  than  that  at  the  anode  (i.e.  Ac  >  Aa). 

3.3.  Validation  of  model  parameters 

The  identification  of  the  impedance  spectra  obtained  by  elec¬ 
trochemical  impedance  spectroscopy  (EIS)  provided  a  set  of 


(a)  Re  Z  (Ohms) 


0.26  0.28  0.3  0.32  0.34  0.36  0.38  0.4  0.42  0.44  0.46 

(c)  Re  Z  (Ohms) 


Fig.  7.  Impedance  spectra  of  the  fuel  cell  versus  its  fractional  model  fitted 
j  =  0.072  A  cm-2,  (d)  j  =  0.074  A  cm-2. 


model  parameters  whose  relevance  must  be  checked.  A  conve¬ 
nient  way  to  do  this  consists  in  comparing  experimental  and 
simulated  Nyquist  plots  of  the  fuel  cell  at  a  given  current  den¬ 
sity  using  the  parameters  identified  in  Section  3.2.  Fig.  7  displays 
the  Nyquist  plots  versus  its  fractional  model  fitted  for  the  current 
densities  of  0.068,  0.070,  0.072  and  0.074  A  cm-2  respectively. 
Then,  it  can  be  observed  that  the  results  of  modeling  are  in  good 
agreement  with  the  measurements. 

4.  Dynamic  response 

In  previous  sections,  we  showed  that  the  insertion  of  a  half¬ 
order  Warburg  impedance  in  the  equivalent  circuit  model  of  the 
fuel  cell  made  it  possible  to  have  a  reliable  and  compact  fre- 
quential  model.  The  fractional  approach  is  then  at  first  place  a 
frequential  approach.  Before  concluding  this  paper,  it  seemed 


0.25  0.3  0.35  0.4  0.45  0.5  0.55 


(b)  Re  Z  (Ohms) 


0.24  0.26  0.28  0.3  0.32  0.34  0.36  0.38  0.4 


(d)  Re  Z  (Ohms) 

measured  (*),  simulated  ( — )  (a)  j  =  0.068  A  cm-2,  (b)  j  =  0.07  A  cm-2,  (c) 
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1.  Fractional  transfer 
function 


2.  Generalized 
differential  equation 


3.  Generalized  system 
of  state  equations 


4.  Solution  of  state 
equations 


Fig.  8.  Steps  for  calculating  the  time  domain  response  of  a  system  containing 
non-integer  derivatives. 


however  essential  to  us  to  transform  this  frequential  model  into 
time  domain,  which  is  an  obligatory  step  for  the  study  of  the 
dynamic  behavior  of  the  fuel  cell. 

The  principal  steps  of  the  generalized  algorithm  for  cal¬ 
culating  the  time  response  of  a  fractional  system  containing 
non-integer  derivatives  are  summarized  in  Fig.  8  [16]. 

Each  step  of  this  algorithm,  which  involves  extensive  calcu¬ 
lations,  is  described  in  Appendix  A.  This  algorithm  has  been 
applied  to  determine  the  time  response  of  the  air-breathing 


v(t) 


Fig.  9.  Simplified  schematic  of  a  PEMFC. 


PEMFC  studied  in  Section  3.  In  order  to  find  the  transfer  func¬ 
tion  of  the  fuel  cell,  the  load  current  i(t)  is  taken  as  input  of 
the  model,  and  the  voltage  across  its  terminals  v(t),  as  its  out¬ 
put.  The  equivalent  circuit  of  the  fuel  cell  is  shown  in  Fig.  9.  The 
reversible  thermodynamic  potential  Et h  (=  1 .23  V)  of  the  fuel  cell 
is  represented  by  a  DC  voltage  source. 

The  fractional  transfer  function  of  an  electrode  is  given  by 


Zk(s)  = 


Vk(s ) 


sC 


eff 

dl 


sC 


eff 

dl 


Rtk(j)  T  Z^k(s) 
+  Rtk(j)  +  Z^k(s) 


where 


AkU)  •  ^/mk 

's/ S  Ci)()k 


(18) 


where  Vk(s)  represents  the  Laplace  transform  of  vk(t ),  and 
co ok  =  1/to k  the  angular  cut-off  frequency.  After  the  application 
of  the  algorithm  described  above,  the  output  voltage  v(t)  of  the 
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Fig.  10.  Simulated  voltage  responses  of  the  fuel  cell  and  the  cathode  for:  (a)  rising  step  input;  (b)  falling  step  input. 
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Table  2 

Values  of  model  parameters  obtained  for  j  =  0.06  A  cm 


Parameters 

Values 

Units 

^ta 

1.94  x  10"8 

Rtc 

3.27  x  10“5 

Aa 

~  0 

£2  rad 1/2 

Ac 

4.858  x  10“6 

Q  rad 1/2 

0.3 

y-Qii 

152 

F 

T 

0.0404 

s 

fuel  cell  is  then  given  by  the  following  equation: 

lit)  =  Et h  -  E-'fZaW/W}  -  £-'{Zc(s)l(s)}  -  Rmi(t)  (19) 

where  I(s)  represents  the  Laplace  transform  of  the  input  current 

at). 

In  order  to  find  the  voltage  response  of  the  fractional  model, 
two  types  of  input  signals  i{t)  were  applied;  a  rising  step  signal 
of  amplitude  1.5A  (Fig.  10a- 1),  and  another  falling  step  signal, 
of  the  same  amplitude  (Fig.  10b- 1).  The  final  results  of  the  sim¬ 
ulation  of  the  voltage  response  of  the  fuel  cell  (Fig.  10,  a-1  and 
b-1)  and  of  cathode  (Fig.  10,  a-2  and  b-2)  are  shown  on  Fig.  10. 

For  the  simulation  of  Fig.  10,  the  values  of  the  parame¬ 
ters  identified  by  EIS  (in  Section  3)  were  linearly  extrapolated 
in  order  to  obtain  the  parameters  corresponding  to  1.5  A  (i.e. 
j  =  0.06  A  cm-2).  These  values  are  given  in  Table  2.  Note  that  the 
transient  response  is  in  good  agreement  with  the  time  response 
of  the  gas  transport  [12]. 

5.  Experimental  validation  of  time  response 

In  this  section,  the  dynamic  response,  derived  from  the  frac¬ 
tional  model  of  PEMFC  in  Section  4,  has  been  experimentally 
validated. 

In  order  to  experimentally  measure  the  dynamic  response  of 
a  commercial  air-breathing  PEMFC  on  the  test  bench,  a  peri¬ 
odic  square  wave  signal  of  the  current  of  amplitude  1.5  A  and 
frequency  0.1  Hz  was  imposed  at  the  fuel  cell  terminals  by  an 
electronic  load.  The  corresponding  variations  in  the  cell  volt¬ 
age  were  then  recorded;  the  results  are  shown  in  Fig.  11.  These 
measurements  were  made  at  the  ambient  temperature  (30  °C). 
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Fig.  11.  Measured  dynamic  response  of  the  air-breathing  PEMFC,  cell  voltage 
( — ),  cell  current  ( — ). 
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Fig.  12.  Comparison  of  measured  and  simulated  dynamic  responses  of  the 
PEMFC,  simulated  ( — ),  measured  ( — ). 

The  results  of  simulation  (Fig.  10)  and  experimental  mea¬ 
surement  (Fig.  11)  are  compared  for  a  period  in  Fig.  12.  It  is 
observable  that  the  modeling  results  are  in  good  agreement  with 
the  measurements.  The  mean  error  over  this  period  is  equal  to 
0.64%,  and  the  maximum  error  is  about  4.13%. 

6.  Conclusions 

A  fine  and  compact  frequential  model  of  a  PEMFC  taking 
account  of  the  diffusion  and  charge-transfer  phenomena  was 
presented  in  this  paper.  This  new  modeling  technique  distin¬ 
guishes  completely  from  the  traditional  approaches  since  it  uses 
non-integer  derivation.  In  particular  it  makes  possible  to  reach 
a  model  of  the  fuel  cell,  reliable  on  a  wide  range  of  frequencies 
with  a  limited  number  of  parameters.  Moreover,  each  parameter 
of  the  model  has  a  physical  significance. 

The  technique  of  electrochemical  impedance  spectroscopy 
was  demonstrated  as  a  powerful  tool  in  order  to  validate  this 
model.  Lastly,  a  general  algorithm  for  calculating  the  time 
response  of  fractional  systems  was  used  for  the  simulation  of  the 
transient  state  of  the  fuel  cell.  In  order  to  validate  the  simulated 
dynamic  behavior  of  the  fuel  cell,  a  commercial  air-breathing 
PEMFC  was  subjected  to  rapid  variations  in  the  load  current 
and  the  output  voltage  was  recorded.  The  results  of  simulation 
were  found  in  good  agreement  with  experimental  results. 

However,  a  further  study  is  necessary  for  the  development 
of  a  more  specific  algorithm  for  the  simulation  of  the  transient 
response.  The  major  work  needed  in  future  consists  in  defining 
more  systematic  methods  for  identification  of  the  parameters 
and  their  validation  particularly  in  the  time  domain.  Moreover, 
it  would  also  be  interesting  to  study  in  detail  the  influence  of 
the  working  conditions  and  geometrical  parameters  of  the  fuel 
cell  (i.e.  geometry  of  electrodes,  the  thickness  and  the  type  of 
membrane  etc.)  on  the  model  parameters. 

Appendix  A.  Calculation  of  time  response  of  the  fuel  cell 

Step  1:  Definition  of  the  transfer  function 

We  can  describe  a  system  in  the  form  of  a  transfer  function 
between  an  output  O  and  an  input  I,  as :F(s)  =  where  O 


1180 


M.  U.  Iftikhar  et  al.  /  Journal  of  Power  Sources  160  (2006)  1170-1182 


Fig.  Al.  Schematic  of  an  electrode. 


Step  3:  Generalized  differential  equation 

The  Eq.  (A. 2)  can  also  be  written  in  the  form: 

Bo  ■  VOi  -  <y0)  +  Bi/2s\/2V(si  -  «o)  +  V(si  “  <A)) 

+  B3/2s\/2V(s\  -  <wo)  =  A0/Ol  -  «o) 

+  Ai/2s\/2I(si  -  a>o)  (A. 3) 

We  know  that  [12]: 


and  I  indicate  the  Laplace  transforms  of  the  functions  o  and  i, 
respectively.  The  objective  is  of  course  to  determine  the  output 
oft)  for  a  given  input  ift). 

Generally,  the  function  F(s)  can  be  written  in  the  form  of  a 
quotient  of  two  polynomials  A(s)  and  B(s),  each  containing  the 
integer  derivatives  of  s  and  one  or  more  implicit  systems  of  frac¬ 
tional  nature  all  of  the  same  angular  cut-off  frequency  coo  =  1/to. 

First  of  all,  we  can  write  an  operational  relation  between  the 
current  I(s)  and  the  voltage  V(s)  for  a  single  electrode  (Fig.  Al): 

y(s)_  f/R/f+f/rf)] 

I(s)  A  +  R^i>  +  zw2w 

where 


Si  ■  G(s  1  -  co o)  =£  {Z)(Q,,[g(f)  •  exp(too)];si) 

where  G  is  the  Laplace  transform  of  the  function  g.  If  we  apply 
the  inverse  Laplace  transform  to  the  Eq.  (A. 3),  the  following 
generalized  differential  equation  is  obtained: 

£0n(Oexp(to0)  +  Bi/2D{l/2]\v(t)  exp (too)] 

+  BiD{l)[v(t)exp(tooo)]  +  #3/2£>(3/2)[f(0  exp(too)] 

=  A0/(Oexp(to0)  +  Ai/2D(1/2)[/(Oexp(too)]  (A.4) 

In  order  to  simplify  the  notations,  we  consider  next  the  variables 
s(t)  and  eft)  such  that: 


A  O')  •  y/cop 

\A  +  COO 


where  j  represents  the  current  density.  In  this  case,  A(s)  and  B(s) 
are  defined  by 

A(s)  =  Rt(J)  +  Z[f(s) 


s(t)  =  vft)  •  exp(too),  eft)  =  ift)  •  exp(too) 


By  dividing  each  term  of  the  Eq.  (A.4)  by  #3/2,  we  can  then 
write 


gV2„(l/2) 

£3/2 


ft)  + 


AL,d )(,)  +  s(3/2)(f)  =  E(t) 

B3/2 

(A.5) 


B(s)  =  1  +  sCdl  •  Rt (j)  +  sCa  •  z\f(s) 


where 


Step  2:  Change  of  the  variable  of  Laplace 

The  second  step  consists  of  carrying  out  the  change  of  the 
variable  of  Laplace,  i.e.,  s\=s  + coo. 

In  effect,  the  transfer  function  (A.l)  becomes  a  rational  frac¬ 
tion  of  polynomials  in  s\  of  only  integer  and  non-integer  orders: 


Aoe(t)  +  Ai/2^2)(t) 
B3/2 


Step  4:  Change  of  variables:  non-integer  derivative  of  the 
output 


Lpi  -  oop)  _  Api  -  cop) 
I(s  1  -  cop)  B(s  1  -  cop) 


Ao  +  A  1/2*54 


1/2 


Bp  +  B\/2S\  +  B\s\  +  B3/2SY 


where 


Aq  = 


Mj)  •  ff 

Cdi 


A 1/2  = 


RtU) 

Cdi 


Bo  =  -  A(j)  •  <ol /2, 


1 

B 1/2  =  — - coo  ■  Rt(j) 

Cdi 


B\  =  A(j )  •  w1/2, 


The  change  of  variables  carried  out  in  this  step  is  the 
key  of  this  algorithm  of  calculating  time  response.  Indeed, 
it  is  the  operation  which  will  enable  us  to  build  the  gener- 
(A.2)  alized  system  of  state  equations,  and  thus  to  solve  the  prob¬ 
lem. 

The  system  (A. 6)  below  clarifies  better  this  change  of  vari¬ 
able: 

s(t)  =  x\ ft),  s^2\t)  =  Da/2)s(t)  =  x^2\t)  =  xz(t), 
S(1\t)  =  Da)s(t)  =  41/2)(0  =  X3(t), 

s(3/2\t)  =  D0/2)sft)  =  41/2)(0  =  x4  (t)  (A.6) 


B3/2  =  Rt  O') 


The  unknown  to  be  determined  is  of  course  x\ . 
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Step  5:  Construction  of  generalized  system  of  state 
equations 

By  using  the  change  of  variables  defined  by  the  system  (A. 6), 
the  Eq.  (A. 5)  can  be  written  in  the  form: 


The  first  term  at  right-hand  side  corresponds  to  what  we  are 
seeking  to  calculate  at  the  instant  tm  (tm  =  m-h),  the  sum  Qj(m-h) 
corresponds  to  the  samples  calculated  at  the  previous  instants. 
Applying  the  expression  (A.  11),  for  each  state  variable  of  the 
system  (A. 9),  we  obtain: 


x^2\t)  =  *4  (t)  =  —a\x\(t)  —  <22*2(0  —  <23*3(0  +  E(t) 

(A. 7) 


with 

Bo  Bin  ,  Bin 

ci\  =  - ,  op  =  — —  and  <22  =  — — 

B3/2  B3/2  B3/2 

From  the  Eqs.  (A. 6)  and  (A.7),  one  can  then  write  the  generalized 
system  of  state  equations  between  the  input  vector  e(t ),  the  state 
vector  x(t)  =  [*i(t),*2(t),*3(t)]T,  and  the  output  vector  s(t): 

X(t/2)(Y)  _  ^  _|_  p  .  s _  Q  .  x (y)  P)  .  e (y) 


where 


0 

1 

0 

"0" 

A  = 

0 

0 

1 

B  = 

0 

— <21 

-<22 

-<23 

1 

0  0 


D  = 


0 


(A. 8) 


Step  6:  Discretization 

In  order  to  solve  the  problem,  we  will  first  of  all  modify  the 
form  of  the  system  of  state  equations;  starting  from  the  relations 
(A.7)  and  (A. 8): 


-1  0 

0 

h -V2  _1 

<2l 

<22  <23  +  h 

- Q\(mh ) 

— 

-Ql(mh) 

—  Qlimh)  +  E{mh) 

x\(mh) 

X2  (mh) 
x3(mh ) 

(A.  12) 


Step  7:  Solution  of  the  system 


By  inversing  the  matrix  M  defined  in  the  relation  (A.  12),  it 
is  then  possible  to  determine  the  value  of  each  state  variable 
at  any  instant  and  particularly  that  in  which  we  are  interested, 
namely  x\.  By  multiplying  this  variable  by  the  term  exp(— t/r0) 
for  every  instant  t ,  we  finally  obtain  the  variable  v(t).  Following 
the  same  steps  for  the  anode  and  cathode,  and  then  using  the 
equation  of  polarization  of  the  fuel  cell,  subtracting  the  ohmic 
drop  in  the  membrane,  cell  voltage  can  be  plotted  with  respect  to 
time. 

uce\\ (0  =  1.23  -  |na(0l  -  \vc(t)\  -  Rm  •  i(t) 

where  pa(0  and  vc(t)  are  the  output  voltages  of  the  anode  and 
cathode  respectively,  determined  using  the  above  algorithm. 


*1  l/2\t)  -  *2(0  =  0,  x(2/2\t)  -  *3(0  =  0, 

-41/2)(0  +  a\x\{t)  +  <22*2(0  +  <23*3(0  =  Bit)  (A.9) 

Then,  we  compute  each  variable  with  a  step  size  of  h.  The  follow¬ 
ing  relation  gives  an  approximation  of  the  non-integer  derivative 
of  a  function  at  a  sampling  instant  tm  =  m-h ,  [12]: 

m 

xf\tm)  =  y2xv3k)  ■  Xj[(m  -  k)-h ]  (A.  10) 

k= 0 


where 


XVj(k)  = 


1 


(-D 


kVj{vj  -  1)  X  •••  X  (Vj  -  k+  1) 


hvik\ 


1 


hvJ 


for  k  >  1 
fork  =  0 


where  vj  indicates  the  order  of  the  derivation  of  the  variable  xj. 
In  the  next  we  will  consider  for  k  >  1 : 

rDjiDj  ~  1)  X  •  •  •  X  (Vj  -  k  +  1) 


covAk)  =  (-  iy 


k\ 


Eq.  (A.  10)  can  then  be  written  in  the  form: 


1 

hvJ 


1 

•  *  /  (fin)  T  - 

J  hvJ 


m 

■  S2a>vj  (kCxj  [( m 

k=  1 


k)-h] 

(A.  11) 
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